Strong correlations via constrained-pairing mean-field theory 
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We present a mean-field approach for accurately describing strong correlations via electron number 
fluctuations and pairings constrained to an active space. Electron number conservation is broken 
and correct only on average but both spin and spatial symmetries are preserved. Optimized natural 
orbitals and occupations are determined by diagonalization of a mean-field Hamiltonian. This 
constrained-pairing mean-field theory (CPMFT) yields a two-particle density matrix ansatz that 
exclusively describe strong correlations. We demonstrate CPMFT accuracy with applications to the 
metal-insulator transition of large hydrogen clusters and molecular dissociation curves. 
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The efBcicnt and accurate description of strong (static 
or nondynamical) correlations in electronic structure the- 
ory has remained an elusive goal. Despite its impor- 
tance in many physical and chemical processes, a first- 
principles, accurate, black-box, and computationally in- 
expensive scheme for static correlation remains unknown. 
We present in this Letter an interesting alternative that 
allows to treat strong correlations in large molecular sys- 
tems and opens a way to treat extended systems too. We 
build upon previous work of Staroverov and Scuseria[l|, 
who showed that a Hartree-Fock-Bogoliubov (HFB)0, 



scheme with an effective pairing interaction of — can 
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accurately describe the strong correlations occurring in 
the 4-electron series and in molecular dissociations done 
over the correct spatial- and spin-symmetry restricted 
HF (RHF) reference. For C = 1 (referred to as IHFB be- 
low), this scheme is equivalent to the functional known 
as Corrected Hartree-Fock (CHF)[3l in density matrix 
functional theory [^. IHFB allows for efficient optimiza- 
tion of natural orbitals (NOs) and occupations via ma- 
trix diagonalization of a mean-field Hamiltonian of small 
dimension (twice the number of orbitals). The HFB 
wavefunction breaks electron number conservation thus 
requiring a chemical potential to adjust the electron 
count to its correct value N. HFB is connected to the 
A^-electron multi configurational antisymmetrized gemi- 
nal power wavefunction [l| (known as pfaffian in recent 
literature^), thus the appearance of static correlation 
should not be unexpected. For most molecules near equi- 
librium, IHFB does not include correlations and yields 
the RHF solution. As the system dissociates, IHFB adds 
strong (left-right) correlations but significantly overcorre- 
lates at large internuclear separations in almost all cases. 
This behavior was considered a serious drawback Q, and 
therefore the model lost attention. In this Letter, how- 
ever, we show that by incorporating simple additional 
constraints, a mixture of HF and HFB methods for in- 
active and active sets of orbitals, respectively, can over- 
come the overcorrelation problem and yield an accurate 
mean- field approximation for strong correlations. Our 
model is here denoted constrained-pairing mean-field the- 
ory (CPMFT). To add important dynamical correlations 



near equilibrium, we propose a straightforward extension 
of the TPSS correlation functional [Sj, which is applica- 
ble to the hydrogen clusters discussed below. Through- 
out this Letter, we will use exclusively the real NO basis, 
{(fii}, where the one-particle density matrix (IPDM) is 
diagonal. 

Theory. We partition the orbitals into core, active, 
and virtual spaces. Only the orbitals in the active space 
are subject to a Bogoliubov transform; the core and vir- 
tual orbitals are treated by HF. In this way, we constrain 
pairing interactions exclusively to the active space. The 
CPMFT wavefunction is therefore a BCS ansatz 0, [|| 
mixed with HF 
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Here |0) indicates the vacuum, a bar on a suffix means 
down spin, and indices c and i run over core and active 
spaces, respectively. Normalization requires xf + yf = 1 
and yf — rii are occupation numbers of the IPDM jij = 
{a\aj). Eq.(IT|) includes iV-electron determinants as well 
as configurations with electron number different from N . 
Double excitations are included in the CPMFT ansatz. 

and 



The pairing matrix K,y = {a\a},) is = 7 ^ 
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= v/n-j(l - n.i). 
_ — Ufj — n and 
or 1 and k,: is ^ 



is diagonal in the NO basis, = Xiyi 
Our formalism is closed-shell, i.e. Ua 
< n < 1. Note that Ki is when n 
when correlation is all static, e.g. at dissociation. More 
details about HFB can be found in Refs.[l|, [3, Q. 

In CPMFT, we introduce additional constraints to di- 
vide orbital spaces into fully occupied (core, n = 1), frac- 
tionally occupied (active), and virtual (n = 0), resulting 
in the energy expression 
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n, - TVa , (2) 



where v runs over the virtual space, Nq and A^a are the 
the number of electrons in the core and active spaces, 
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respectively, and the different /x's are the chemical po- 
tentials (Lagrange multipliers) constraining the electron 
numbers in each space. The IHFB energy above is 



(3) 



where hij, Jij 



{ipi(pj\(p,(pj), and K^j = {ip,ipj\ipjip,) 
arc one-electron, two-electron Coulomb and exchange in- 
tegrals, respectively. The two-particle density matrix 
(2PDM) of this ansatz is 



ij 



1 



{liklji ~ laljk - KijKki) , 



(4) 



which includes particle-particle and hole-hole but no 
particle-hole correlation terms. As shown below, this 
simple ansatz describes strong correlations very accu- 
rately. We conjecture that the cumulant term (k Cg) k 
with = 7 — 7^) offers a compelling definition of what 
we intuitively understand as static or strong correlation. 
Note that the negative sign of the last term in Eq.Q 
comes from the attractive pairing potential. 

The mathematical solution to the CPMFT constrained 
problem is to set fic = —00 and fiy = 00. The physi- 
cal meaning of this is that /ic (/^v) pulls down (pushes 
up) the energies of core (virtual) orbitals so that these 
levels do not mix with the active orbitals. In our cal- 
culations, we set the chemical potential for the inactive 
spaces /iv = ~MCi a-nd gradually increase fiy smoothly 
to convergence. 

As a benchmark example, we discuss first the dissocia- 
tion of H2 , a paradigmatic example of strong correlation. 
The CPMFT(2,2) (two electrons in two active orbitals) 
wavefunction for such a system is given by 
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(|0) + l'^<Ti,<^<Ti3) + I'Pal^'P-y'J 



(5) 



Note that the first and last determinants of Eq. ^ con- 
tains and 4 electrons, respectively. However, the ex- 
pectation value of the electron number operator on this 
CPMFT wavefunction is 2. The energy associated with 
Eq.® is 
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Jaiscrl^ Kai^a'l^, (6) 



and given that at dissociation, Jo-iao-*^ = ^o-iao-*^, this 
expression correctly reproduces the energy of two isolated 
hydrogen atoms. Eq.® also yields the exact IPDM at 
dissociation. 

Results. We have implemented CPMFT in the Gaus- 
sian suite of programs [9j. As another paradigmatic ex- 
ample of strong correlation, we present results for the 
dissociation potential of the N2 molecule in Fig. [1] This 
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FIG. 1: Potential energy curves of N2 calculated with a Gaus- 
sian 6-311++G** atomic orbital basis set. 



is a challenging case because it requires up to six-electron 
excitations to yield the correct curve. CPMFT(6,6) gives 
a very accurate description of this dissociation at low 
(mean- field) computational cost. It accomplishes this 
feat while preserving spatial and spin symmetries. Unre- 
stricted HF (UHF) improves over RHF, but the density 
is symmetry-broken with three a electrons localizing on 
one N atom, and three [3 electrons on the other. Com- 
plete Active Space Self- Consistent Field CASSCF(6,6) 
calculations [lO|, which include all excitations for 6 elec- 
trons within 6 active orbitals, also yield a correct poten- 
tial curve for N2 including a part of dynamic and static 
correlations near equilibrium and only strong correlations 
towards dissociation. Note how IHFB overcorrelates over 
the entire curve yielding very poor results. 

While CPMFT significantly improves over HF and 
IHFB in diatomics, we would like to present here its ap- 
plication to large systems that are currently completely 
out of the reach of other approaches for accurately de- 
scribing electron correlations. Accordingly, we will inves- 
tigate below large hydrogen networks. Before doing so, 
it should be evident from Fig. [T] that CPMFT lacks dy- 
namical correlation, especially at equilibrium, where the 
total correlation energy is underestimated. Density func- 
tional theory (DFT) seems a natural candidate for adding 
this effect. However, it is known 1 1 . [T^ that semilocal 



DFT correlation functionals face substantial challenges 
when applied to symmetry adapted densities as those re- 
sulting from CPMFT (as opposed to symmetry broken 
unrestricted models). Alternatives to circumvent this is- 



sue have been proposed in the literature [ill Il3l |. In or- 



der to deal with this problem and with large hydrogen 
clusters in mind, we here propose a simple solution that 
uses K as a "detector" of static correlation. We introduce 
a correction factor to DFT correlation functionals using 
gi^l- {2Ki)'^ and define g{r) = J2i 5i"i|v'j(r)P, where 
m is a parameter controlling the rate of change of gi. We 
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FIG. 2: Potential energy curves for H2 dissociation calculated 
with a Gaussian cc-pV5Z basis set. 



thus write 



,(r).9(r)dr, 



(7) 



where ec(r) is the DFT correlation energy density. We 
have modified the TPSS correlation functional to include 
g{r) and define E^'^^^^. 

In Fig. [21 we plot the potential energy curve of II2. For 
Rh-h 00, correlation in this system is all static while 
near equilibrium, there is substantial dynamical correla- 
tion. Here, £:^tpss jg evaluated with the CPMFT density 
and NOs, hence our calculations are not self-consistent, 
and we shall denote them as CPMFT(kTPSSc). The pa- 
rameter m is chosen as 10 throughout this Letter based 
on a fit of CPMFT(kTPSSc) to the the exact full config- 
uration interaction (FCI) [lO| dissociation curve of H2 . 

We now discuss our results on hydrogen networks, 
namely, a one dimensional (ID) H50 chain and a three 
dimensional (3D) 6x6x6 hydrogen cube, with an inter- 
nuclear distance Rh-h between adjacent atoms. In the 
limit of Rh-h — > 00, these model systems dissociate to 50 
and 216 isolated hydrogen atoms, respectively. As Rh-h 
increases, they display a metal-insulator transition. Both 
systems are par adig matic models for strongly correlated 
Mott insulators [15j and are not treatable by CASSCF 
due to the staggering number of active configurations, 
10^^ and 10^^'^ for the ID and 3D cases, respectively. 
In the linear c ase, the density matrix renormalization 
group (DMRG)[lJ] approach, a very accurate wavefunc- 
tion method, is available [l3|- However, for the 3D case, 
DMRG is not (yet) applicable. As shown below, accurate 
CPMFT of these two model systems only requires diago- 
nalizing Hamiltonians of moderate dimension, 100 x 100 
and 432 x 432, respectively. 

Fig. [3] presents potential energy curves for each system 
computed with several methods and an ST0-6G basis. 
This choice of minimal basis is made to allow compari- 
son with results in the literature tul. All electrons and 
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FIG. 3: ( Top) Potential energy curves for the symmetric disso- 
ciation of a H50 chain. DMRG results are taken from Ref . [l5| . 
(Bottom) Same as Top, but for a 6 x 6 x 6 hydrogen cube. 



orbitals are explicitly treated in our calculations. The ac- 
tive spaces of CPMFT are (50,50) and (216,216) for the 
ID and 3D systems, respectively. This means that for 
these clusters Nc = and Na is equal to the size of the 
active space in Eq.(l2]). In both cases, RHF and second- 
order M0ller-Plesset perturbation theory fMP2)(loj yield 
unreliable curves, while coupled cluster singles, doubles, 
and perturbative triples (CCSD(T))[l^ has convergence 
difficulties at long internuclear distances and hence is not 
included in the plot for the 3D case. The RHF method 
misses a considerable amount of strong correlation. MP2 
may both miss or exaggerate non-dynamical correlation, 
as goes forth from Fig. 3. CCSD(T), despite its sin- 
gle reference character, is quite good at covering non- 
dynamical correlation, thanks to its inclusion of infinite- 
order terms in a balanced way. If there are problems 
for long R values, they are not visible in Fig. 3 due to 
the mentioned convergence problems. For the ID system 
around equilibrium, the DMRG energy is very similar 
to that of CCSD(T). While CCSD(T) breaks down past 
the equilibrium bond length, DMRG is able to predict 
accurately the potential curve toward dissociation. 
On the other hand, CPMFT and CPMFT(kTPSSc) 
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TABLE I: Correlation energy (a.u.) of a H50 chain at Ro 



Basis set 


MP2 


CCSD(T) 


CPMFT CPMFT(fi:TPSSc) 


STO-GG" 


-0.545 55 


-0.765 47 


-0.211 12 


-1.384 31 


cc-pVDZ 


-0.976 45 


-1.213 41 


-0.203 42 


-1.426 27 


cc-pVTZ 


-1.123 26 


-1.347 9'' 


-0.205 46 


-1.408 85 


cc-pVQZ 


-1.162 74 


-1.374' 


-0.205 81 


-1.410 06 



"The DMRG correlation energy with STO-6G basis is -0.772 67 
(Ref.lli). 

''Estimated value by extrapolaion of a H„ chain to n = 50. 
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FIG. 4: Decay of off-diagonal density matrix term (712) for 
two hydrogen atoms at diagonal vertices in a 4 x 4 x 4 hydrogen 
cube. 



yield correct dissociation curves (within the basis set 
used) not only for the ID system but also for the 3D case. 
Only CPMFT properly describes the metal-insulator 
transition of the 3D hydrogen network. Although at first 
glance CPMFT(kTPSSc) appears to overcorrelate the ID 
chain near equilibrium, this is partially a basis set effect: 
in the top panel of Fig. [3] CCSD(T) is essentially under- 
estimating correlation effects due to the small basis used 
in the calculation (a necessity to make the CCSD(T) cal- 
culation affordable). 

To further investigate the accuracy of CPMFT, wc 
have carried out calculations of the ID hydrogen chain 
near equilibrium (Rc = 1.8 bohr) with very large cc-pVnZ 
basis sets, where n = D,T, Q. We compare our correla- 
tion energies with CCSD(T), which yields results close 
to DMRG (Fig. [m with the ST0-6G basis and thus we 
consider it very accurate for the ID chain near Re. For 
the cc-pVTZ and cc-pVQZ bases, we extrapolated the 
CCSD(T) correlation energy of H50 from shorter chains, 
as it was not feasible to compute CCSD(T) for H50 within 
reasonable time. Results are presented in TABLE H] 
While convergence of the CCSD(T) correlation energy 
with basis set size is slow, the correlation energies of both 
CPMFT and CPMFT(/tTPSS) converge reasonably well. 
Static correlation is expected to be fairly basis set inde- 
pendent, a property reproduced by CPMFT. 



In order to visualize the metal-insulator transition ex- 
plicitly, we have plotted the off-diagonal density matrix 
term 712 between two hydrogen atoms at diagonal ver- 
tices in a 4 X 4 X 4 cube (Fig. S]). Both RHF and MP2 off- 
diagonal terms remain substantially above zero as Rh-h 
increases, implying that the electrons are delocalized. 
The CPMFT off-diagonal 712, on the other hand, rapidly 
decays toward zero. Evidently, only CPMFT reveals the 
gradual metal-insulator transition in this hydrogen cube. 

Finally, we would like to mention some current limita- 
tions and prospects of our CPMFT model. When applied 
to hetero-nuclear systems, CPMFT docs not reproduce 
the exact dissociation energy although there is a very 
substantial improvement over RHF and IHFB. For ex- 
ample, the CPMFT dissociation energy error relative to 
the correct CASSCF hmit is 0.2 eV for BH with a 6- 
311+-f G(3df,3pd) Gaussian basis set. Ways to address 
and correct dissociations in hetero-nuclear systems to 
non-degenerate orbitals are currently being investigated 
and will be discussed in future publications. A more 
elaborate scheme for dealing with dynamical correlations 
for symmetry adapted densities is also being developed. 
Most importantly, we would like to stress our conjecture 
that the particular CPMFT form of the 2PDM in Eq.® 
is most adept for describing strong correlations. 
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